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by Leslie L. Scalzott and Carl F. Lorenzo 
Lewis Research Center 

SUMMARY 

A stability criterion has been developed for use with digital simulations or graphical 
data. The criterion used is a simple one, which basically requires ultimate state bound- 
edness of the time function representing the system output. Techniques are established 
so that a transient of finite duration may be used as a practical indication of stability. 

A stability test utilizing the criterion has been implemented into a general digital 
program which will allow the determination of system stability directly by computer. 

This program and criterion was demonstrated for three applications representing: (1) a 
linear system, (2) a continuously nonlinear system, and (3) a discontinuously nonlinear 
system. 


INTRODUCTION 

Currently, there is an increased interest in simulation of physical systems using 
digital computers. This is especially true of highly nonlinear systems and systems of 
great complexity both linear and nonlinear. Of particular interest in such systems is 
the response of the system to either a change in operating point or response to some 
form of input disturbance. Generally it is desirable to form a stability map showing the 
effect of two or more system parameters on stability, where stability is some desirable 
property defined by the investigator. 

Presently, to obtain a stability map of a simulated system, a great many transients 
must be simulated in some organized method, all these transients must be examined to 
determine whether or not the response is stable by some criterion, and either these re- 
sults must be mapped directly or further refined with more transient responses. Clearly 
it would be desirable to eliminate reading all these transient responses out of the digital 
computer to form the stability map. Indeed, if a computer program could be evolved 



which would indicate the stability of the system, the user would be saved by the burden of 
applying the stability criteria to the transients and the time he must wait between running 
separate transients. Thus a stability map could be found directly, without the investiga- 
tor in the loop. 

A substantial background of material exists in the literature in the general area of 
stability. References 1 and 2 are examples of papers that deal with the general problem. 
Most of the work had as its objective the analytical determination of stability for systems 
of differential equations. Notable, of course, are the classical efforts of Liapunov, 
Poincare, and Lagrange. Current efforts in the English literature (much of the effort 
has been by Russian investigators) are those of Bellman (ref. 3) and LaSalle and Lefshetz 
(ref. 4). 

The above efforts, as stated, are directed toward analytical solutions of the stability 
problem. To the authors' knowledge, no work has been done in implementing a stability 
criterion that can be applied directly to the numerical solutions of systems of differential 
equations. Direct implementation of these classical criteria is not practical, since it 
would require an unperturbed solution and the establishment of the classical e - 5 rela- 
tion (necessitating a large number of e trials). (Symbols are defined in appendix A. ) 
Reference 2 presents a discussion of this classical criteria. Instead, the present study 
requires the solutions to have certain desirable properties which approximate the char- 
acter of those classical criteria. These required properties are then formalized into a 
working stability criterion. Such a stability criterion and its implementation into a digi- 
tal program are the subjects of this report. 

The results of this study should also be useful in the areas of digital adaptive control 
systems, data analysis, and as a working criterion for use with experimental studies. 


GENERAL SCHEME FOR STABILITY ANALYSIS 


The general scheme for determining the stability of a time function generated by 
digital simulation can be understood by reference to figure 1. There are three basic 
elements under consideration: (1) the simulation proper, which generates an output time 



Figure 1. - Block diagram for computer stability analysis. 
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function in serial form, (2) a routine, which examines the time function and determines 
its stability, and (3) a control logic routine, which would change the system parameters 
to achieve the neutral stability points. 

Since a practical approximation to stability is being developed, it is clear that a 
transient of infinite duration cannot be examined. Hence, the stability analysis routine 
must be capable of making a judgment as to system stability based on a finite observation 
time. The overall scheme, therefore, will be to assign an observation time over which 
the time function will be analyzed and determined to be either stable, unstable, or neu- 
trally stable. If the observation time is not long enough, the time function will be condi- 
tionally stable or conditionally unstable. If a conditional situation develops in the search 
for a stability map, the control logic also has the function of progressively increasing the 
simulation time by a time increment until an unconditional stability or instability occurs. 
The criteria on which the above conditions are based will be defined formally in the next 
section. 

This report basically deals with the stability analysis block. Throughout the studies 
presented herein, relatively simple control logic is used. The part of the control logic 
which determines the new system parameters can be complex. This problem is strongly 
analagous to the problem of determining roots of equations using digital computer tech- 
niques, for which an excellent background of material is available (e. g. , ref. 5). 

The overall problem can be reduced to the following: Given a function of time which 
represents the output of a system to some bounded input disturbance, determine the sys- 
tem stability in some finite observation time. 

Clearly the choice of a value for observation time is critical in that the user must 
recognize for his system the dynamic elements with the longest response time and choose 
the observation time appropriately. 

In the preceding discussion only a single system output has been examined. If the 
possibility of separate parts of a simulation independently being stable or unstable exists, 
then all these outputs must be monitored for stability. 


STABILITY CRITERION 

The previous section discussed the general problem of digital determination of the 
stability of a time function. It was premised on the existence of a technique to discern 
stability using only the perturbed transient responses. The stability analysis block 
(fig. 1) should be capable of handling any type of time function in order to be generally 
useful; that is, time functions should include responses which are continuous or discon- 
tinuous, oscillatory or nonoscillatory, and periodic or aperiodic. Furthermore, the sys- 
tem responses of both linear and nonlinear systems are of interest. The analysis should 
be flexible in order to have general applicability and should also be easily implemented 
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on the digital computer for use with digital simulations and should use computing time 
economically. 

Fundamentally the criterion is an approximation to asymptotic stability; that is, the 
requirement will be ultimate state boundedness of the output time response. 

Stability criterion fo r infinite observation time. - A system is con- 
sidered to be stable if, after the completion of a bounded perturbation 
on the system, its motion never exceeds some e distance from the 
bias after some time tp that is, 

|Bias - f(t)| < e for t > t^ 

where 

Bias = lim — f f(t)dt 
T-oo T T) 

If there does not exist a time tj, such that the motion never exceeds 
the prescribed e distance from the bias, the system is called un- 
stable. 

Physically, the stability definition requires that, after the system transients (due to 
either a perturbation or change in operating point) die out, the system response remains 
within a prescribed bound about the new operating point. It further considers simple 
divergences and growing envelopes to be unstable. 

If the stability of linear constant coefficient systems based on the roots of the char- 
acteristic equation is considered, the above definition of stability is compatible with lin- 
ear stability for all points with the possible exception of those on the imaginary axis; 
that is, a neutrally stable linear system may oscillate forever at some amplitude greater 
than the prescribed e bound and, based on the given definition, would be unstable. This 
will not be a problem, however, since, in the study of such systems, the neutral stability 
point (linear stability) will generally be found by approaching from the positive and nega- 
tive damping sides and will, therefore, yield the same stability map. 

The following is a practical criterion based on the preceding infinite observation 
time definition. 

Stability criterion for finite observation time. - A system is con- 
sidered to be unconditionally stable if, after the completion of a 
bounded perturbation on the system, the motion never exceeds some e 
distance from the bias after some time tj but before completion of 

some observation time t . It is further required that the slope of the 

1 / 

stability curve be negative. The system is conditionally stable if, 
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under the preceding conditions, the slope of the stability curve is posi- 
tive. If there does not exist a time t., < t such that the motion does 

1 c 

not exceed the prescribed e distance and if the slope of the stability 
curve is positive the system is unconditionally unstable. If, under 
these conditions the stability curve has a negative slope, the system 
is conditionally unstable. The stability curve is defined as the locus 
of | Bias - Extrema | . 

This criterion is also stated in table I. The difference between the finite and infinite 
observation time cases is the additional slope requirement for the stability curve. This 
slope requirement is used as an indication of the future convergence of the envelope of 
the time function. 


TABLE I. - THEORETICAL STABILITY CRITERIA FOR 
FINITE OBSERVATION TIME 


Functions 

Condition 

Stability curve 
boundedness requirement 

Slope 

requirement 

Oscillatory 

Stable, unconditionally 

|Bias - Extrema| < e 

Slope < 0 


Stable, conditionally 

|Bias - Extrema| < e 

Slope > 0 


Unstable, unconditionally 

|Bias - Extrema| > e 

Slope > 0 


Unstable, conditionally 

|Bias - Extrema| > e 

Slope < 0 


Neutrally stable 

|Bias - Extrema| < e 

Slope = 0 


Neutrally unstable 

|Bias - Extrema] > e 

Slope = 0 

Nonoscillatory 

Stable, unconditionally 
Unstable, unconditionally 

Bias converges 
Bias does not converge 



From a practical point of view, this definition still presents some difficulties. Con- 
sider, for example, a case of a time function with beats. While the function may be 
bounded, it may be beating very slowly, and give an apparent positive slope for the sta- 
bility curve. This difficulty, however, can be handled by programming techniques which 
consider the extrema of the envelope instead of the extrema of the time function. This 
will be described in detail in the next section. 

In essence this finite observation time criterion, modified by the above technique, is 
the stability analysis which has been implemented in this study. Further checks and fea- 
tures that have been built into the program are discussed in the next section. 

This criterion is different from the classical criteria of Liapunov, Poincare, and 
Lagrange as discussed in reference 2, since a different goal is desired. The classical 
cases generally compare the perturbed and unperturbed responses. Also this is done 
analytically instead of numerically. In the present situation the stability is determined 


5 



I I 


III II II II 


I ii ii mi 


without necessarily having the mathematical form of the response and without having the 
unperturbed motion of the system . 


DIGITAL PROGRAM 

A digital computer program was written based on the criterion presented in the pre- 
vious section. This program is called DIGSTA (digital stability analysis). A flow dia- 
gram of the computer program and a complete FORTRAN IV listing are given in appen- 
dix B. All the checks and definitions applied to the system under investigation are illus- 
trated in this appendix. DIGSTA is written in FORTRAN IV computer language which can 
be readily translated into other languages. A simplified flow diagram of DIGSTA is pre- 
sented in figure 2 to assist in understanding the general logic of the stability routine. 

This figure also illustrates the general scheme used for the determination of stability. 
Here, the following logic pattern is used. 

The time function generated by the system response is used as an input f(t). A 
check is made to see if the bias value has been found. If the bias is not found, a time 
check is made. If t < t , the investigation continues with a new f(t) point. If the bias 
does not converge for t > t , the response cannot be bounded and is therefore considered 

V 

unstable (see proof in appendix C). This assumes that t is chosen to be a large enough 
finite time to be representative of an infinite time for the system being studied. If the 


From MAIN 
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bias is found, a check is made on boundedness for some prescribed number of extrema. 

If the boundedness condition is satisfied, the system is stable and the slope of the stabil- 
ity curve is then checked. If the slope is negative, the system is unconditionally stable; 
if it is positive, it is conditionally stable; and if it is zero, it is called neutrally stable. 

If the system is not bounded for N extrema in one observation time, a slope check is 
made again. Herein, positive slopes indicate unconditional instability and negative slopes 
conditional instability. 



Figure 3. - Typical transient. 

There are several terms used in DIGSTA which should be defined or mentioned here. 
Figure 3 shows a basic curve which could be operated on by DIGSTA. This figure iden- 
tifies names referred to in this report and in the routine. 

Before a boundedness check can be applied, the bias of the system must be estab- 
lished. The bias level has been attained when 


_d_ 

dT 


1 _ 

T 



= 0 


In the program the running average (RAV) is defined as 


n 

RAV s _J_ V f.(t)At n = — 
nAt^j 1 At 
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The bias is the converged value of the running average and, therefore, ideally requires 

— (RAV) = 0 
dT 

Two possible numerical requirements for convergence are 

iRAVj -RAVjJ < |RAV i _ 1 |AtC 

and 

|RAV i - RAV U1 [ < AtC 

The first condition is satisfactory when the RAV is a large number since it requires the 
running average to converge to within some percentage of the required value. When the 
bias is near zero, however, RAV is approaching zero. Therefore, the condition re- 
quires a smaller and smaller percentage and convergence is not achieved. The second 
condition is satisfactory around zero; however, large values of RAV could require too 
small a percentage error. The program, therefore, uses condition 1 for 

RAV > 1 


and condition 2 for 


RAV < 1 


As mentioned previously, the slope of the stability curve is used in defining the 
system as being conditionally, unconditionally, or neutrally stable or unstable. In 
DIGSTA only the sign of the slope is of interest and this sign is obtained by an analysis 
of the signs of each side of the envelope of the function. The sign of one side is defined 
as the arithmetic sign of (EXTR. - EXTR. „) and is called SLOPE1. The sign of the 

other side is defined as the arithmetic sign of (EXTR^ j - 
EXTR i 3 ) and is called SLOPE2. If sign of (EXTR^ - 
EXTR^j) is positive (negative) then SLOPE 1 is the upper 
envelope of the curve and SLOPE2 is the lower envelope 
and conversely if (EXTRj - EXTR^j) is negative. The 
sign of the slope of the stability curve then is given by 
SGN(MM SGN(EXTRj - EXTR^j)) where the dummy vari- 
able MM is determined from table n. 


TABLE H. - DUMMY VARIABLE 


AS FUNCTION OF SLOPE 


Slope 2 

Slope 1 

+ 

- 

0 

+ 

±9.9 

-1.0 

-1.0 

- 

+1.0 

±9.9 

+ 1.0 

0 

+1.0 

-1.0 

0 
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Now, a -1. 0 value of MM SGN(EXTR^ - EXTR^_j) indicates that the envelope is 
converging and a 1. 0 value indicates divergence. Then this product together with the 
boundedness requirement gives a stability indication. 

The +9. 9 value denotes that the signs of the envelopes are contradictory (i. e. , both 
positive or both negative) and is a warning to the investigator. Herein, +9. 9 indicates 
an unbounded function and -9. 9 indicates bounded. 

DIGSTA performs a beat and high frequency oscillation investigations on the func- 
tion of interest. This is done by examining the envelope of each side of the time func- 
tion for changes in the arithmetic sign of the slope. Changes in sign are an indication 
of those phenomena. When this occurs, DIGSTA refers to the extrema of the extrema 
(BEXTR) in lieu of EXTR. This method is illustrated graphically in figure 4. 

Since DIGSTA is a routine used in conjunction with a main program, certain con- 
stants and variables must be made common to both programs. These terms are 
IPRINT, K1K2, L, M, N, DELTAT, EPSLON, T, TPREV, Y, YP, TCRIT, and 
TINCR. They are defined in appendix A and values are recommended where applicable. 

Certain discretion must be exercised in defining some of the constants. In particu- 
lar, N must be chosen so as to be sufficiently large to account for peculiarities en- 
countered when investigating a system that is beating or has high-frequency oscillations 
superimposed; that is, N must be larger than the number of consecutive, similarly 
signed slopes of odd-numbered intervening extrema (see fig. 5). 

In addition to being the number of successive stability curve points required to be 
within an e-band of the bias for a system to be defined as being stable, N is also used 
as 
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(1) The number of successive slopes required to define a beat or high frequency 

oscillations 

(2) The number of successive slopes required to define a system as no longer 

beating or having high frequency oscillations (2N) 

If the user desires to alter any of these criteria, he can readily do so by adjusting the 
pertinent FORTRAN "IF statement" in DIGSTA. 

Similarly, the user must have an insight into his problem so he can select reasonable 
TCRIT, TINCR, DELTAT, and EPSLON values. Adequate time must be allowed to per- 
mit the program to solve for the bias and then apply the system stability investigation. 

In most applications the settling time can generally be gaged by examination of those 
system components with the longest characteristic time. For first order systems the 
time constant t is the controlling factor, and for the second order, the parameter 

is the controlling factor. A choice of TCRIT of 10 or more times the greatest 
characteristic time should be suitable for most "linear like” systems. Smaller e val- 
ues require larger TCRIT values. 

For an unstable configuration this amount of time should be adequate to establish the 
instability. The critical factor here is that the user at least allows adequate time 
(TCRIT) to solve for the bias for the stable situation. 

It is the authors' feeling that the BIAS as presently defined is satisfactory for most 
purposes although somewhat slow to converge. If, after TCRIT has elapsed, the system 
has been identified as being conditionally stable or unstable, the user may deem it appro- 
priate to increase the observation time. Therefore, TINCR may be assigned a nonzero 
value. Also, DELTAT (At) must be chosen to be sufficiently small to reasonably define 
the extrema. Using information concerning the system under investigation, the user de- 
fines an EPSLON value in the main program which establishes the maximum acceptable 
amplitude for system stability. 

One merit of DIGSTA is the small amount of storage space required (113 storage 
locations). Also, as a typical example, the relay controller simulation and its logic sec- 
tion required 64 storage locations. This total is small relative to that available on the 
larger computers. This advantage will prove particularly useful when DIGSTA is coupled 
with a main program that requires a large number of storage locations. Also, DIGSTA 
can be easily implemented to the program of interest since either equations or data can 
be operated on by the stability program . 

To handle physical data containing noise, some care is necessary in the selection of 
the N number as discussed previously. Another asset of DIGSTA is the small additional 
time requirement. An exponentially decaying cosine function was operated on by DIGSTA 
and was found to be stable after 0. 19 minute of computer time and 18. 101 seconds of 
simulation time (DELTAT = 0.001 sec). This transient was then run without DIGSTA for 
18. 101 seconds of simulation time. The computer time required for this case was 
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0. 09 minute. This example should not be interpreted to mean that the use of DIGSTA 
doubles the time required. The transient routine could be quite complex, and it could 
have a large computer time requirement. This example merely demonstrates that 18 101 
data points were analyzed by DIGSTA in 0. 10 minute. This may be used to estimate the 
time required by DIGSTA. 


Control Logic 

Some control logic is needed to adjust the system parameter being varied to achieve 
the stability locus. In the applications that follow an elementary technique was used to 
iterate to the neutrally stable condition. Two values were chosen to be the bounds of the 
parameter. One of these was assumed to be stable and the other unstable. An initial 
parameter value between the two limits was chosen. If the system was then found to be 
stable, a constant value was added to the initial guess until a parameter value was found 
that rendered an unstable system (and, conversely, for the unstable system). Then it 
was necessary to iterate between these last two values. The method used was that of 
choosing a value midway between these extremities and thereby establishing new bounds. 
This process was repeated until the two limits were an assigned percentage apart. Then, 
the neutrally stable point was defined as the stable limit. This technique can be modified 
or replaced in order to increase efficiency. One alternative for linear-like systems 
would be to use the slope of the envelope of the last trial as a basis for selecting a new 
trial point. 

For the computer used, there is a restart capability. This means that, in the event 

- 38 

that the adjustable parameter is so chosen that there is an underflow (X < 10 ) or an 

90 ~~ 

overflow (X > 10°°) or a division by zero (computer arithmetic errors), the control logic 
section can then adjust the parameter and the investigation can be continued. DIGSTA 
has been written to allow reentering the stability subroutine after such a computer oper- 
ation stoppage without a loss of continuity. 


APPLICATION OF DIGITAL PROGRAM 

The program developed in the preceding section has been applied to a series of ele- 
mentary functions to verify program operation and has also been applied to three cases 
of interest to demonstrate use of the program. 

The elementary functions were sines and cosines with exponential decay and growth 
functions as multipliers. These were studied with and without bias and combinations of 
functions were studied to simulate a beating situation. All these functions gave proper 
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stability indications with the criterion as previously defined. 

The three practical cases include study of (1) the stability of a linear system, (2) a 
nonlinear system with a continuous type (mild) nonlinearity, and (3) a system with a dis- 
continuous type (hard) nonlinearity. 


Linear System Stability 

A stability study was made of the dynamics of a simple monopropellant rocket engine 
system with a linear injector. It consists of an injector of the form Wj = K(P,p - P^), a 
pure dead time representing propellant vaporization time and a first order lag represen- 
tation of the chamber fill dynamics. The equations defining this system are as follows: 

dP r 

T — + p c = K > b 
dt u 1 D 

W b (t) = w,(t - a) 

*i = K 2 < p t - P C> = K 2 AP 

These equations are expressed in block diagram form in figure 6, using Laplace notation. 



Figure 6. - Block diagram for rocket engine stability study. 

Numerical approximations to the above equations were used to generate the time function 
required for the stability analysis. 

The following nominal conditions were used: 


9 fi 

Chamber pressure differential, P^,, psi; N/m 682; 4.70x10° 

Weight flow, w, lb/sec; kg/sec 463; 210 

Chamber gain constant, Kj, sec/in. ; sec/cm 1.4721; 0.2282 

Delay time, a, sec 0.002 

Time constant, r, sec 0. 001 
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Figure 7. - Digital implementation for rocket engine chugging problem. 
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The critical time t c for this analysis was chosen as 1 second. Also N = 6 and a 
At of 0. 00001 second were used. A flow diagram which shows how the simulation, sta- 
bility subroutine and control logic were combined is presented in figure 7. For 
e = 1500 psi (10. 34x10 N/m^), the program was asked to determine the injector AP to 
within 1 psi (6.89x10 N/m ) required for stability for a given dead time. This, in turn, 
was repeated for the number of dead time values required to form the graph of figure 8. 
The disturbance used to excite the system was a step in the tank pressure P T from 0 to 
1130 psi (0 to 7. 79X10 6 N/m 2 ) at time zero. 

The area above the curve in the figure indicates injector pressure drops for which 
the system is stable, and that below the curve, unstable. It will be noted that stability 
as determined from the roots of the characteristic equation (see ref. 6 for technique) fell 
on the same line within readability limits. 

To get the 12 stability points shown in this figure required the simulation of 156 
transients. The e value of 1500 psi (10. 34X10 6 N/m 2 ) was used in the stability criterion 
since the initial large transient would have taken significant time to arrive within a 
smaller e band. This is permissible since linear system stability is being studied. 


Nonlinear System Stability (Continuous Nonlinearity) 

In order to study a mild nonlinearity, the rocket system simulation discussed in the 
above section was modified by assuming a nonlinear injector characteristic; that is, 
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Wi - k 2 y^p-p - p c 

Also, (P,p - P^-,) was limited to zero so that the back flows were not possible. The nomi- 
nal conditions are the same as the above case. For these conditions, e values of 10 and 
100 psi (6. 89xl0 4 to 68. 9xi0 4 N/m 2 ) were studied. The stability of this system in re- 
sponse to the same input as studied previously is shown in figure 9. 



Normally, to handle this type of nonlinearity, the method of small perturbations 
would be applied to linearize the system equations. These, in turn, are analyzed based 
on the roots of the characteristic equation. This analysis has been done and the results 
of such a linearization are also plotted in figure 9. The difference is an indication of the 
effect of the nonlinearity on system stability. The effect of the square law injector is to 
distort the weight flow wave shape. This distortion shifts the average-value weight flow 
during dynamic oscillations, which, in turn, increases the pressure drop required to 
stabilize the system. 

The magnitude of the bias shift depends on wave amplitude. Thus, the magnitude of 
the difference between the stability curves is dependent on the size or harshness of the 
input disturbance. 

In this study the smaller values entailed longer simulation times than would be en- 
countered in the linear case. 


Nonlinear System Stability (Discontinuous Nonlinearity) 

For this application, a stability analysis was made for a system in which a relay is 
used to control the position of an inertial torsional servo. The block diagram for this 
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Figure 10. - Block diagram for position control torsional servo, 
system is shown in figure 10. The defining equations for this system are 


T 




-<<P 

2 

~-<(p <* 

2 ~ ~ 2 

, a 

(p <-- 
2 

<P = e set - e 

The relay controller is being studied, to demonstrate the use of DIGSTA for limit cycle 
forms of operation. The nominal conditions for this analysis were 


Damping constant, K n , ft-lb-sec; m-kg-sec 0.0001; 0.138x10"^ 

Moment of inertia, I, ft-lb-sec; m-kg-sec 1.000; 0.138X10 

Limit of angular position, a, rad 0.01745 

Desired Q value, 0 se j., rad 0 

In the DIGSTA subroutine the following values were used: 

Critical time, t , sec 50 

Number of extrema, N 5 


+1 


T R = < 


V.- 
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In this particular study, some difficulty was encountered in arriving at a valid simu- 
lation. As a result, since the problem is piecewise linear, the differential equations 
governing the behavior of the output over any time interval were solved analytically. 

This yielded an algebraic equation which in turn was solved to generate the system out- 
puts. This equation was solved in a serial form, which formed an adequate simulation 
input for the digital subroutine. The relay deadband was ±0. 0087265 radian (±1/2°). 

The disturbance was an initial angular velocity of 0. 005 radian per second. 

The results of this study are shown in figure 11. Here, for a given value of damping 
a Kj (within 0. 001) was established which would indicate stability for some value of e. 

In this case, the e value corresponds to the acceptable limit cycle amplitude. Also in- 
dicated in this study are the effects of variations of servodamping K^. 

The study involving limit cycles presented some small difficulty since long simula- 
tion times were required to arrive at final stability indications. For K D = 1x10"^ and 
e = 0. 00878 radian, 44. 7 seconds of transient were required. This is due to the slow 
convergence or divergence of the envelope inherent in this problem. 

The above are typical applications in linear and nonlinear controls and dynamics 
systems. The program could also be used with a real system as an indication of stabil- 
ity. This is possible since the time function generated from the computer simulations 
are in serial form, and the subroutine would not see any difference. Some problems in 
control logic would exist. However, it might not even be necessary to have a control 
logic unit if conditional stability indications were acceptable for a given application. 
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This ability to examine real systems opens the avenue to possible applications in 
digital adaptive control systems. 


SUMMARY OF RESULTS 

The most important result of this work is that a stability criterion has been evolved 
which will allow a stability analysis of systems for bounded, finite duration inputs. The 
criterion is general in that it will function properly for practical time functions, and it 
is flexible to the extent that the user can pick his boundedness band and special index 
numbers. The techniques evolved in the present study can also be used as a basis for 
different or more sophisticated criteria. The problem of working with a finite observa- 
tion time has been satisfactorily handled by use of the stability curve, that is, the infor- 
mation in the envelope of the time function. 

The stability criterion has been implemented for use with the digital computer and 
requires modest additional computing time to determine the stability bounds directly, 
furthermore, a FORTRAN program implementation of the criterion together with a 
simple control logic has been presented. 

The program has been demonstrated for three applications which cover a linear, 
continuously nonlinear, and discontinuously nonlinear dynamic systems. For these sys- 
tems stability maps have been generated directly from the computer without the need of 
examining separate system transients. 

The program allows the reduction of turn-around time in the use of the digital com- 
puter for mapping stability results, and allows a uniform basis for stability analysis. 

The control logic used in the applications presented has been simple, it is likely that 
a significant economy in computing time can be gained by implementing more sophisti- 
cated techniques, that is, techniques analogous to those used for root finding of equations. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, June 29, 1967, 

180-31-01-05-22. 
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APPENDIX A 


SYMBOLS 


Engineering Symbols: 
a 8 limit, rad 

C convergence constant 

2 2 

I moment of inertia, lb -ft -sec ; kg-m-sec 

K gain constant 

Kj chamber gain 

Kg injector gain 

P pressure, psi; N/sqm 

AP pressure differential, psi; N/sq m 

s Laplace variable 

T torque 

t time 

At time increment 

t c critical time 

w weight flow, lb/sec; kg/sec 

6 angular position 

0 ge t desired 8 value 

cp position error 

e allowable oscillation amplitude 

0 delay time 

r time constant 

£ damping constant 

natural frequency, rad/ sec 

Subscripts: 
b burning 

C chamber 

D damping 

1 injector 

i current 


R required 

T tank 

1 initial 

FORTRAN Variables Common to Both the Main Program and the Stability Routine: 

DELTAT time increment; defined in MAIN 

EPSLON allowable oscillation magnitude; defined in MAIN 

IPRINT when IPRINT = 0, there is a printout every stable or unstable transient; if 
IPRINT * 0, there is a printout only for a neutral point (slope = 0); de- 
fined in MAIN 

K1K2 identifies the system as being stable, unstable, conditionally, uncondition- 
ally or neutrally stable or unstable. Utilized by subroutine change logic; 

calculated in DIGSTA 

L equal to 1 when f(t) is not equal to constant for first time; set in MAIN 

M equal to 0 until DIGSTA is entered; set in MAIN 

N number of successive stability curve points required to be within e-band of 

bias; defined in MAIN 

T current system running time; developed in MAIN 

TCRIT allowable time to identify the behavior of the system; defined in MAIN 

TINCR allowable increment to increase TCRIT; defined in MAIN 

TPREV previous T; defined in MAIN 

Y variable being tested in DIGSTA; developed in MAIN f(t) 

YP previous Y; defined in MAIN 

Other FORTRAN Variables Used in the Stability Program: 

BIAS converged value of the running average 

EXTR extrema 

KOUNT1 number of successive points (BIAS) identical within e^-band 

KOUNT2 number of successive stability curve points within e-band of BIAS 
RAV running average 
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diagram (fig. 12) and a listing of the program are presented. 



Figure 12. - Flow diagram of computer program. 
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SLOPE OF SLOPE of EXTR of 


Conditional! 



Unconditionally 


/ ls\ 

/ time \v 



System is 

conditionally 

stable 


If X<0 
XEXTR / 


System is 


No stability 
curve points 
could be found 


System is 
neutrally 
stable or 
unstable 


' If X<0 

^SLOPE/^ 

System is 
conditionally 
stable or unstable 


System is 
unconditionally 
stable or unstable 


Figure 12. - Concluded. 




r> a o o o o 


SUBROUTINE DIGS! A ( IPR INT , K 1K2 , l ,M, N, DELTAT , EPSLGN, T , TPRE V , Y , YP , 

1 TCR IT, T I NCR ) 

DELTAT, EPSLON, TCR I T , T I NCR , IPRINT, AND N MUST BE INITIALIZED IN THE 
MAIN PROGRAM 

T, TP, Y, YP, L, AND M MUST BE CALCULATED IN THE MAIN PROGRAM 

LOGICAL A,B,C,D,E,F 
DIMENSION XA(1I) 

IF (M .GT. 0) GO TO 20 

READ ( 5,7000) ( XA ( II ) , 1 1 = 1 , 1 1 ) 

WRITE (6,7000) ( X A < 1 1 ) , 1 1 = 1 , 1 1 ) 

20 D = .TRUE • 

IF (L .GT. 1) GO TO 40 

C 

C A TRUE - SEARCHING FOR MINIMIM 

C B FALSE - TO RESET KOUNTl = 0 

C C = TRUE - BIAS HAS BEEN FOUND 

C D FALSE ~ EXTREMUM HAS BEEN FOUND 

C E TRUE - NO BEAT PHENOMENA OR HIGH FREQUENCY OSCIL. HAVE OCCURRED 

C F TRUE - NO SECOND ORDER BEAT PHENOMENA OR HIGH FREQUENCY OSC. 

C IPRINT = 0 - INTERMEDIATE PRINT OUTS ELIMINATED 

C J NO. OF EXTREMA 

C JJ NO. OF EXTREMA OF EXTREMA 

C K NO. OF S.C. POINTS 

C KOUNTl = NO. OF SUCCESSIVE PTS. (BIAS) IDENTICAL WITHIN El BAND 

C K0UNT2 = NO. OF SUCCESSIVE S.C. PTS. WITHIN E BAND OF BIAS 

C K0UNT3 = NO. OF SUCCESSIVE SLOPES INDICATING A BEAT OR HIGH FREQ. OSC. 

C KQUNT4 = NO. OF SUCCESSIVE SLOPES INDICATING NO BEAT OR HIGH FREQ. OSC. 

C K0UNT5 = NO. OF SUCCESSIVE SLOPES INDICATING A SECOND ORDER BEAT OR 

C HIGH FREQ. OSC. 

C K0UNT6 = NO. OF SUCCESSIVE SLOPES INDICATING NO SECOND ORDER BEAT OR 

C HIGH FREQ. OSC. 

C K1K2 = K1 + K2 , INDICATES NATURE OF STABILITY OF THE SYSTEM 

C L =1 WHEN Y IS NOT ZERO FOR THE FIRST TIME 

C M =0 WHEN THE SUBROUTINE IS ENTERED FOR THE FIRST TIME 

C N NO. OF SUCCESSIVE S.C. PTS. REQUIRED TO BE WITHIN E BAND OF BIAS 

C BEXTR = EXTREMA OF EXTREMA 

C EXTR = EXTREMA 

C RAV = RUNNING AVERAGE OF DATA POINTS FOR THE SOLUTION OF THE BIAS 

C RX VALUE OF RUNNING AVERAGE WHEN KOUNTl = 0 

C 

C XEXTR = STABILITY CURVE (S.C.) POINT (LOCUS OF EXTREMA) 

C INITIAL CONDITIONS FOR RUNNING AVERAGE AND STABILITY CURVE CHECKS 

C 

BIAS = l.OE+35 

BEXTR = 0.0 

BEXTP = 0.0 

BEXTPP = 0.0 

BEXTP3 = 0.0 

C6XTR = 0.0 

CEXTP = 0.0 

CEXTPP = 0.0 

CEXTP3 = 0.0 

EXTR = 0.0 

6XTRP = 0.0 

EXTRPP = 0.0 

RAV = 0.0 

RX = 0.0 
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SLOPE 1 = 0.0 

SL0PE2 = 0.0 

TB = 0.0 

TBP = 0.0 

TB PP = 0.0 

TC = 0.0 

TCP = 0.0 

TC P P = 0.0 

TT = 0.0 

TTP = 0.0 

T TP P = 0.0 

TR A V = 0.0 

X =0.0 

XEXTR = 0.0 

XRAV = 0.0 

XY = 0.0 

Y 1 Y 

Z = 0.0 

IK =0 

J =0 

JJ =0 

JB =0 

JBP = 0 

JS =0 

JSP = 0 

K =0 

K 1 =0 

KIK2 = 0 

K0UNT1 = 0 

K0LJNT2 = 0 

K0UNT3 = 0 

K0UNT4 = 0 

K0UNT5 = 0 

K0UMT6 = 0 

B = .FALSE. 

C = .FALSE. 

E = .TRUE. 

F = .TRUE. 

C 

C RUNNING AVERAGE CHECK 

C 

40 RPP = RP 

50 RP = RAV 

IF <T .EO. 0.0) GO TO fO 

RAV = (RP*TPREV/DELTAT + Y)*DELTAT/T 

IF ( ABS ( RAV - RX) .LE. XY ) GO TO 70 

C = .FALSE. 

BIAS = l.OE+35 

TBIAS = 0.0 

K0UNT1 = 0 


70 

IF 

(L 

- 2) 

1150, 

80, 

80 

Y 2 


= Y 




IF 

( Y2 

.GT. Yl) 

GO TO 

100 


A 


= .TRUE. 




GO 

TO 

1150 



100 

A 


= .FALSE. 




GO 

TO 

1150 




C STABILITY CURVE CHECK 


110 
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non 


130 


no 

IF 

(Y - YP) 

120, 

650 

120 

IF 

(A) 

GO 

TO 

650 


A 

= .TRUE. 





GO 

TO 150 




130 

IF 

(A) 

GO 

TO 

140 


GO 

TO 650 




140 

A 

= .FALSE. 




150 

J 

= J + 1 





EXTRP3 = EXTRPP 
EXTRPP = EXTRP 
EXTRP = EXTR 
EXTR = YP 

BEAT OR HIGH FREQUENCY OSCILLATION CHECK 

EXEX = EXTR - EXTRPP 
JSPP = JSP 
J SP ~ j s 

IF ( A B S ( EXEX) .GT. 1 . 0E-07*EXTR ) GO TO 170 

JS = JSPP 

GO TO 180 

170 JS = SIGISM 1.0, EXEX) 

180 IF (J .LT. 4) 

IF (J/2-2 .NE. J) 

IF (E) 

GO TO 650 

210 IF (JS .EQ. JSPP) 

K0UNT4 = 0 

K0UNT3 = K0UNT3 + 1 
IF ( K0UNT3 • LT • N) 

IF (IK .LT. 2 ) 

GO TO 300 

250 IF (JS .EQ. JSPP) 

K0UNT4 = 0 

IF (E) 

GO TO 310 

270 IF (E) 

K0UNT4 = K0UNT4 + 1 
IF (K0UNT4 .LT. 2*N) 

E = .TRUE. 

K0UNT2 = 0 

K0UNT3 = 0 

IK =0 

GO TO 550 

290 K0UNT3 = K0UNT3 + 1 
IK = IK + 1 

IF (K0UNT3 .LT. N) 

300 E = .FALSE. 

JSPP = JS 
XYZ5 = EXTRPP 

310 IF (JS) 

320 XYZ5P = AMAX1 (XYZ5tEXTKP) 

XYZ5 = AH INI ( EXTR , EXTRP, EXTRPP , EXTRP 3 , X YZ 5 ) 


330 

IF (JSPP) 

XYZ5P = AMINKXYZ5, EXTRP) 

XYZ5 = AMAX1 ( EXTR, EXTRP, EXTRPP, EXTRP3,XYZ5) 

650, 

650, 

340 

340 

IF (JSPP) 

REXTP3 = BEXTPP 

340, 

650, 

650 


BEXTPP = BEXTP 
BEXTP = BEXTR 


GO TO 550 
GO TO 250 
GO TO 210 

GO TO 550 

GO TO 550 
J = J+l 

GO TO 270 

GO TO 290 

GO TO 550 

GO TO 310 


GO TO 550 


320, 650, 330 
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BEXTR = XYZ5P 

K0UNT3 = K0UNT3 + 1 

JJ = JJ + 1 

BX8X = BEXTR - BEXTPP 

JBPP = JBP 

JBP = JB 

IF UBS(BXBX) .GT. 1 .0E-07*BEXTR ) 

JB = JBPP 

GO TO 370 

360 JB = SIGN! 1 .0,BXBX) 

370 IF (JJ .LT. 4) 

IF (JJ/2*2 .NE. JJ) 

IF (F) 

GO TO 650 

400 IF (JB .EG. JBPP) 

K0UNT6 - 0 

K0UNT5 = K0UNT5 + 1 
IF ( KQUNT5 .LT. N) 

IF ( K I .LT. 2) 

GO TO 490 

440 IF (JB .EG. JBPP) 

K0UNT6 = 0 

IF (F) 

GO TO 500 

460 IF (F) 

K0UNT6 = K0UNT6 + 1 
IF (K0UNT6 .LT. 2*N) 

F = .TRUE. 

K0UNT2 = 0 

K0UNT5 = 0 

K I =0 

GO TO 540 

480 K0UNT5 = K0UNT5 + I 
K I = K I + I 

IF (K0UNT5 .LT. N ) 

490 F = .FALSE. 

JBPP = JB 
XYZ6 = BEXTPP 

500 IF (JB) 

510 XYZ6P = AMAX1 ( XYZ6,BEXTP ) 

XYZ6 = AM INI ( BEXTR tBEXTPt BEXTPP t BEXTP3 t XYZ6 ) 
IF (JBPP) 

520 XYZ6P = AMIN1(XYZ6,BEXTP) 

XYZ6 = AMAXI ( BEXTR tBEXTPt BEXTPP *BEXTP3*XYZ6) 

IF (JBPP) 

530 CEXTP3 = C EXT PP 


C E XT PP 

= 

CEXTP 

CEXTP 

= 

CEXTR 

CEXTR 

= 

XYZ6P 

K0UNT5 

= 

K0UNT5 

XYZ1 


CEXTR 

XYZ2 

= 

CEXTP 

XYZ 3 

= 

C EXT PP 

X YZ4 

= 

CEXTP3 

TCPPP 

= 

TCPP 

TCPP 

= 

rep 

TCP 

= 

TC 

TC 


TPREV 

TPPP 

= 

TCPPP 

TP P 

= 

TCPP 


GO TO 360 

GO TO 540 
GO TO 440 
GO TO 400 

GO TO 540 

GO TO 540 
JJ = JJ+I 

GO TO 460 

GO TO 480 

GO TO 540 

GO TO 500 


GO TO 540 

510, 650, 520 

650, 650, 530 

530, 650, 650 
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TP 

= TCP 


TRAV 

= TC 


GO TO 

560 

540 

XYZ1 

= BEXTR 


XYZ2 

= BEXTP 


XYZ3 

= BEXTPP 


XYZ4 

= BEXTP3 


TBPPP 

= TBPP 


TBPP 

= TBP 


TBP 

= TB 


TB 

= TPREV 


TPPP 

= TBPPP 


TPP 

= TBPP 


TP 

= TBP 


TRAV 

= TB 


GO TO 

560 

550 

XYZ1 

= EXTR 


XYZ2 

= EXTRP 


XYZ 3 

= EXTRPP 


XYZ4 

= EXTRP3 


TTPPP 

= TTPP 


T TPP 

= TTP 


TTP 

= TT 


TT 

= TPREV 


TPPP 

= TTPPP 


TPP 

= TTPP 


TP 

= TTP 


TRAV 

= TT 

560 

D 

= .false. 


650 IF (C) 

IF (B) 

RX = RAV 

KOUNT1 = 0 

B = .TRUE, 

GO TO 690 

660 IF (ABS(RAV) .GT. 1.0) 

XY = DELTAT-.01 

GO TO 670 

665 XY = .001*ABS(RX) 

670 IF (ABStRAV - RX > .GT. XY ) 

K0UNT1 = K0UNT1 + 1 
IF (K0UNT1 - N) 

680 B = .FALSE. 

690 IF (T .LT. TCRIT) 

WRITE (6,3000) T , J , RP , R A V , RPP , R X , EXTR , EX TRP 
K 1 =2 

IF (F) 

JZ = JB 

JZP = JBP 

GO TO 720 
710 IF (E) 

•JZ = JS 

JZP = JSP 

720 WRITE (6,6000) X YZ3 , XYZ2 , XY Z 1 
IF ( ABS( JZ) - ABS( JZP) ) 

730 IF ( AB S ( R P - RAV) .GE. ABS(RPP - RP ) ) 

740 K2 =3 

GO TO 1140 
750 K2 =6 


GO TO 780 
GO TO 660 


GO TO 665 

GO TO 680 
690, 760, 

GO TO 1150 

GO TO 710 

GO TO 730 

740, 750, 

GO TO 750 


760 


750 
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GO TO 1140 




760 

C = .TRUE* 

TBIAS = T 




780 

IF (ABS(RAV) .IE. . 01 *ABS ( XYZ 1 ) ) 
BIAS = RAV 

RAV 


0.0 


IF (T .GE. TCRIT) 

GO 

TO 

830 


IF (D> 

XEXTR = ABSIXYZ1 - BIAS) 

GO 

TO 

1150 


IF (XEXTR .GT. 1 .OOOODEPSLON ) 
K0UNT2 = K0UNT2 + 1 

GO 

TO 

810 


IF (K0UNT2 - N) 

1150, 

880, 

810 

K0UNT2 = 0 

GO TO 1150 




830 

IF (K0UNT2 .LE. 0) 
K2 =3 

GO TO 875 

GO 

TO 

850 

850 

IF (XEXTR .LE. 0.0) 
K 1 =2 

GO TO 890 

GO 

TO 

870 

870 

K2 = 6 




875 

WRITE ( 6 , 4000 ) J , KOUNT 2 , T , B I AS , T B I AS , TR AV ♦ TP 
K 1 = 1 

GO TO 1140 




880 

K 1 = 1 




890 

SLOPE 1 = XYZ1 - XYZ3 

IF (ABS(SLOPEl) .LE. 1.0E-5*ABS(XYZ1) ) 

GO 

TO 

900 


IF (SLOPED 

910, 

900, 

900 

Ml =0 

GO TO 930 




910 

Ml = 1 

GO TO 930 




920 

Ml - 2 




930 

SL0PE2 = XYZ2 - X YZ4 

IF ( ABS ( SL0PE2 ) .LE. 1 .0E-5*ABS ( XYZ2 ) ) 

GO 

TO 

940 


IF ( SLOPE 2 ) 

950, 

940, 

940 

M2 =0 

GO TO 970 




950 

M2 - 3 

GO TO 970 




960 

M2 =6 




970 

MM = Ml + M2 + 1 

MN = S IGN( 1 .0 tXYZl-XYZ2 ) 

GO TO ( 980,990, 1000,1000*1010, 1000,990,990, 1010) 

, MM 



980 

SLOPE = 0.0 

GO TO 1020 




990 

SLOPE = -MN 
GO TO 1030 




1000 

SLOPE « +MN 
GO TO 1030 




1010 

XKK = Kl*(l - Kl) + 1 

SLOPE = SIGN( 1 . 0 ,XKK ) *9.9 
GO TO 1060 




1020 

K2 =9 

K 3 = 10 

K4 =11 

GO TO 1110 




1030 

I F ( K 1 .GT. 1) 

GO 

TO 

1100 


IF (SLOPE .GT. 0.0) 

GO 

TO 

1060 

1050 

K 2 =6 

K3 =7 





880 


920 


960 
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K4 =8 

GO TO 1110 
1060 K2 = 3 

K 3 =4 

K4 =5 

IF ( K 1 • GT • 1 ) GO TO 1110 

IF (T - TKAV .LT * T I NCR ) GO TO 1110 

WRITE (6*5000) T * T I NCR 
1090 K2 =6 

GO TO 1140 

1100 IF (SLOPE) 1060* 1060, 1050 

1110 IF (IPRINT .GT. 0) GO TO 1130 

IF ( K2 .LE. 8) GO TO 1140 

1130 WRITE (6*2000) X A ( K 1 ) * SLOPE * B I AS * T * X A ( K2 ) * X A ( K3 ) , X A ( K4 ) * XA ( K 1 ) * 
1XEXTR 

1140 K1K2 = K1 + K2 
1150 RETURN 

2000 FORMAT ( ////20X1 OH SYSTEM IS A2 * 20HST ABLE AND SLOPE = F6.1* 10H. 
1BIAS = G1 5 • 8 * 5X7HT I ME = G1 5 . 8 //45X 10HSYST EM IS 3A6 * A2 * 6HST ABL E / / 

2 50X8HXEXTK = G15.8) 

3000 FORMAT ( // //20X31HBIAS COULD NOT BE FOUND WITHIN G15.8,5H SEC.*5X* 

1 I 5 * 2 1 H EXTREMA WERE FOUND ./ 1 5X5HRP = G 1 5 . 8 * 5X6HR A V = G15.8,5X 
26HRPP = G 1 5 # 8 * 5X5HRX = G 1 5 . 8 /40X7HEXTR = G 1 5 . 8 * 5X8HEXTR P = G15.8// 
325X40HTHE SYSTEM IS CONSIDERED TO BE UNSTABLE.) 

4000 FORMAT (20X*I5,13H EXTREMA AND I5*46H STABILITY CURVE POINTS COULD 
1 BE FOUND WITHIN G1 5 . 8 , 4HSEC . / 1 5X * 7HB I AS = G15.8*10H* TBIAS = 

2G1 5 • 8 * 1 OH * TEXTR = G15.8*19H, PREVIOUS TEXTR = G15.8//8X116HIF THE 

3 PERIOD IS RELATIVELY LONG* INCREASE TCRIT. ALSO* THE BIAS CONVERG 
4ENCE REQUIREMENT MAY BE MADE MORE STRINGENT.) 

5000 FORMAT ( // / / 3 1 X7HT I ME = G15.8*38H SEC. NO EXTREMA HAVE BEEN FOUND 
1 FOR G1 5 • 8 * 5H SEC .//31X69HASSUME THE SYSTEM IS OVERDAMPED AND THER 
2 E FORE UNCONDITIONALLY STABLE.////) 

6000 FORMAT ( 20X93HTHE SYSTEM IS BEATING OR HIGH FREQUENCY OSCILLATIONS 
1 OCCURRED. THE EXTREMA OF THE EXTREMA ARE/ 20X8HEXTPP = G15.8,10X 
27HEXTP = G15.8*10X7HEXTR = G15.8) 

7000 FORMAT (1X2A2*9A6) 

END 

$ DATA 

UN CONDITIONALLY UNCONDITIONALLY NEUTRALLY 
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APPENDIX C 


RUNNING AVERAGE CONVERGENCE 

Theorem : If f(t) is a bounded real function defined on the interval [0, °°) and piece- 
wise continuous, then the running average 



dt 


must converge in the sense that 


lim — (RAV) = 0 
T— 00 dt 


Proof : The running average is differentiable; hence, 


_d_ 

dT 



f(t) ^ _ f(T) _ 1 
T T T 



m 

T 


dt 


(Bl) 


Consider the first term on the right side of equation (Bl). To establish lim of this 

X-^oo 

function, let C > 0 be given. Then, since f(t) is bounded over the interval, there 
exists a number M > 0 such that 


Now, for T > 0, 


I f (T) | < M 


for 0 < T < 00 


|f(T)| 

T 


f(T) < M 
T “ T 


For any C > 0 there exists a number N 


such that, for T > N > 0, 


— < C 
T ~ 
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or 


f(T) 

T 


< C 


Thus, 


lira = 0 (B2) 

T— oo T 

Considering the second term on the right side of equation (Bl), let a number C > 0 be 
given. Again since f(t) is bounded over the interval, there exists a number M such that 

-M < f (t) < M 


It follows therefore that 


Thus, 




M dt = M 
T 


is bounded, and 



for 0 < T < 00 


Now, for T > 0, 
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Again for any C > 0 there exists a number N such that for all T > N > 0 


M< c 

T - 


or 


Thus 



Now, from equations (Bl), (B2), and (B3) it follows that 


lim 



(B3) 


which establishes the theorem. 

Corollary: It logically follows then that if the running average 


RAV 



does not converge in the sense that 
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lim — (RAV) = 0 
T— °° dT 


and, if f(t) is a piecewise continuous real function defined on the interval [0, °°), then 
f(t) cannot be bounded over the interval. 

This, then, establishes the concept that, if given a sufficient amount of time (infinite, 
if necessary) and if the running average does not converge to give the bias, then f(t) is 
unbounded and, therefore, unstable. In DIGSTA this is approximated by requiring con- 
vergence within TCRIT. (Since TCRIT can be picked to look like infinity for most prac- 
tical systems. ) The program further requires that 

— RAV = 0 
dT 


for N consecutive times. 

With slightly more difficulty, the following theorem can also be shown. 

Theorem: If f(t), f(t), . . . , f n (t) are piecewise continuous, bounded real functions 
defined on the interval [0, °°), then the running average 


RAV e 



T 


must converge in the sense that 


lim 

T— °o 


d - (RAV) = 0 
dT n 


for n = 1, 2, . . . 


This is, of course, a stronger convergence requirement. However, if this kind of con- 
vergence is not achieved, then either f(t) or one of its derivatives is not bounded (there- 
fore unstable). This is a stronger stability statement. The convergence requirement of 
the program more closely resembles this requirement (because of the n consecutive 
dRAV/dt = 0 checks). 
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